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Many networks are embedded in physical space and often interact with it. This interaction can be 
exemplified through constraints exerted on network topology, or through interactions of processes defined 
on a network with those that are linked to the space that the network is embedded within, leading to complex 
dynamics. Here we discuss an example of such an interaction in which a signaling agent is actively 
transported through the network edges and, at the same time, spreads passively through space due to 
diffusion. We show that these two processes cooperate or compete depending on the network topology 
leading to complex dynamics. 

Biological, social, and physical networks can be natural or man-made. Understanding the relationship 
between their structure and their dynamics is of great scientific interest 1 " 3 . Some of these networks can 
be studied in isolation because they interact very weakly with the space or the environment that they are 
embedded within. A citation network is one example of such a network - the physical location of the authors does 
not affect the citation of their article. Other networks, however, are strongly affected by the physical space they sit 
in. These are usually referred to as spatial or geographical 4 networks. They are subject to constraints stemming 
from the interaction with the space. This interaction can take different forms, the most common one being 
constraints, dependent on the metric of space, imposed on the network topology. This could be a cost function 
which creates longer connections or a distance dependent probability of connecting to other network elements 5-7 . 
Another type of spatial interaction can be defined when the network topology is independent of spatial (envir- 
onmental) constraints, but the dynamics within the network involve other processes that are linked to or defined 
within the embedding space. In these instances transport mechanisms on the network interact with those defined 
on the environment. Other examples come to mind: neuronal or astrocytic networks that are coupled through 
neurites/gap junctions/synapses, but also secrete chemicals into the extracellular space, or epidemiological net- 
works where the disease spreads through personal contact and through diffusion of pathogens through the air 8 . 

In this paper we study an example of such an interaction. We define a chemical agent which spreads actively 
through the network but is also secreted into the "environment" where it undergoes passive spread governed by 
diffusion. Thus the transport is composed of an active component which allows for (nearly) instantaneous spread 
of the agent through the network, with the macroscopic pattern of its spread being dependent upon the network 
connectivity, and the passive spatial spread outside the network, through "environmental" diffusion. The inter- 
action of these two transport processes, as we will show below, may be highly nonlinear and result in complex 
spatio-temporal dynamics of the network nodes. 

Results 

We define a network composed of simple excitable integrate-and-fire elements (see Methods and Fig. 1 for a 
detailed description). When the level (or concentration) of an agent in the element exceeds a threshold, the 
element rapidly releases the agent from its internal stores (Fig. lb) into the network and the surrounding physical 
space (Fig. la). The transport on the network is an active process - the signal spreads instantaneously to other 
nodes connected with the activated node through standard diffusive coupling. At the same time, the signal 
spreads passively via a spatial diffusion process, the nature of which is determined by the spatial extent and 
properties of the physical space, Fig. la. For simplicity, we consider that both these mechanisms only have positive 
signs - i.e they both spread excitation. This system loosely corresponds for example to astrocytic networks in the 
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Figure 1 | Description of the model, (a) Cartoon highlighting bimodal 
interactions from the active element (left) to the charging element (right). 
The processes underlying this interaction are 1) the active pathway 
through the network, defined by connectivity parameters (connectivity 
radius R, rewiring parameter p) and diffusive network coupling constant /?, 
and, 2) the passive environmental diffusion, where AT amount of the agent 
is emitted into the embedded space, diffuses subsequently with speed D* 
over to the charging element, which absorbs it with proportionality 
constant y. (b) The time dependence of the release of the agent from 
internal nodal stores occurring after the node reaches the concentration 
threshold, (c) An example of the time evolution of agent levels within a 
single node; the small short spikes represent instantaneous changes in 
agent concentration due to noise, after the threshold (dashed line) is 
reached the node releases large amounts of agent into the network and its 
immediate environment. 



brain. Astrocytes, in simple terms, release calcium stored in the 
endoplasmic reticulum into the cytosol, which increases free calcium 
levels in the cell 9 . This calcium spreads through gap junctions to 
other astrocytes, but at the same time activates complex chemical 
cascades responsible for the release of glio -transmitters into the 
extracellular space. Glio -transmitters in turn can stimulate other cells 
(neurons or astrocytes) in the spatial vicinity 1011 . 

While the specific dynamics of network elements can take differ- 
ent forms, we show that, depending on the connectivity of the net- 
work, the active network transport and the passive environmental 
transport sometimes cooperate in spreading the excitation through- 
out the network and at other times compete, effectively inhibiting 
propagation through the system. 

We investigate the spatio-temporal pattern formation within the 
network as a function of the network topology, the strength of net- 
work connections, and the feedback amplitude between the envir- 
onment and the network nodes. It is well established that if the 
magnitude of either network coupling or environmental feedback 
amplitude is large enough then that process alone can drive the 
formation of sustained large-scale spatio-temporal patterns 1213 . 



However, we observe that for intermediate coupling values of both 
processes a strong dependence on network topology emerges. For 
local network topologies the two transport mechanisms cooperate, 
forming globally propagating waves, while for random topologies the 
processes compete, impeding signal propagation throughout the sys- 
tem. This effect is shown on Fig. 2, which shows the mean nodal 
frequency and mean pairwise phase coherence (MPC) between net- 
work nodes as a function of the magnitude of the feedback between 
the environment and the nodes (y), and efficacy of network connec- 
tivity (/?) for four different network topologies. When the network is 
effectively disconnected (/? — > 0) and y < 0.1 1 the mean frequency of 
nodal activity is low as it is driven only by external noise. When a 
certain amplitude (y = 0.11) of nodal interaction with the envir- 
onmental agent is reached the network enters a high activity global 
activation state (Fig. 2). However, if the amplitude of network coup- 
ling is increased for local network topologies (p = 0) the transition to 
sustained global activation of the network takes place for lower values 
of y (Figs. 2a, 2c) - the two transport processes cooperate. In contrast, 
for global network topologies (p = 1) sustained network activation is 
significantly delayed as y increases (Figs. 2a, 2c) - the two transport 
processes compete, impeding each other. Figures 2c and 2d summar- 
ize the difference in network activation for the two topologies as a 
function of both the magnitude of network transport and the envir- 
onmental feedback. The increase of network coupling can play 
reverse roles in affecting the global network activation depending 
on the network topology. 

For local network topologies the cooperation of the two processes 
is due to the fact that both processes support local signal propagation 
from a recently activated node to its neighbors. Active network trans- 
port is faster (it is instantaneous) than environmental diffusion. 
However, on short spatial scales, the excitation stemming from the 
two processes can partially coincide providing supra-threshold 
excitation and leading to activation of new nodes in the proximity 
of the active site. As this process continues away from the originating 
site it forms large scale propagating waves within the network. The 
formation of local patterns that propagate through the network 
depends on the interplay between the magnitude of the diffusion 
constant, D, and the signal decay constant a (see Eqn. (1) in 
Methods section). The transition to global activation patterns as a 
function of the magnitude of network coupling is depicted in Fig. 3d, 
where we show co-activation patterns between network nodes giving 
global and local changes of mean phase coherence and cross-correla- 
tion in network activation sites. For weak network coupling (Fig. 3, 
location i), diffusion itself can not sustain large enough areas of local 
activation for the signal to spread globally through the network. 
Instead, random local domains are formed and then dissipate 
quickly. This leads to relatively low correlation between network 
nodes (Fig. 3b). The mean size of the formed domains can be inferred 
from the distance dependence of cross correlations between the 
nodal activation patterns (Fig. 3c). Characteristic network activity 
plots show the evolving size of synchronous excitations (Fig. 3d, left 
row). Histograms of pairwise cross-correlation values are Gaussian 
due to the stochasticity of the activated locations (Fig. 3d, row i right). 

For intermediate levels of network coupling (locations ii and iii in 
Fig. 3), the domains spread through the network forming first irregu- 
lar large scale activation patterns (Fig. 3, location ii), and then organ- 
izing into propagating waves (Fig. 3; location iii). Large irregular 
areas of activation lead to slower spatial decay of the cross-correla- 
tions (Fig. 3b), and a significant skew toward positive cross-correla- 
tion values on the pairwise histogram (Fig. 3d, row ii right). 
Formation of regular propagating waves through the network 
(Fig. 3, location iii), is characterized by strong correlation of nodal 
activity in the direction normal to the wave propagation, and anti- 
correlation in the direction of wave propagation (Fig. 3b, d; location 
iii). These regular waves result from the front of activation closely 
following nodes which have recently recovered from the previous 
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Figure 2 | Changes in mean nodal frequency (a) and mean phase coherence (b) as a function of relative strength (amplitude) of both transport processes 
(x-axis - strength of interaction of the nodes with the physical space y, y-axis - efficacy of network connectivity /?) and network topology. Every location 
represents results averaged over 9 simulation runs, (c) Changes of nodal mean frequency for local (blue) and random (green) network topologies as a 
function of passive transport amplitude (y); = 0.04. (d) Changes of mean frequency as a function of the active transport amplitude (/?) for local (blue) 
and random (green) network topologies; y = 0.12. 



cycle of activity. The match between the nodal refractory time and 
the speed of the wave propagation also explains the observed non- 
monotonic behavior of the mean nodal frequency (Fig. 2a). Further 
strengthening of the coupling of either process increases the speed of 
network activation, which in turn leads the excitation to self-extin- 
guish due to the nodal refractory time. The activation needs then to 
be randomly restarted which leads to a reduction of the mean nodal 
frequency. The non-monotonicity resulting from the size of the 
environment does not change the overall observation that global 
network activation occurs at lower coupling strengths for local net- 
work connectivities as opposed to higher coupling values for random 
ones. 

Finally, when the network coupling is high the whole network 
activates quickly (Fig. 3d; location iv). The location initiating the 
burst and the pattern of spread are random. Furthermore, anti- cor- 
related islands may form if a given location was recently activated 
because of the noise and is in its refractory state. 

For random network topologies (rewiring probability p = 1) 
increases in network connectivity strength have the opposite effect. 
Even for relatively strong interactions with the environment, the 
network is unable to support global activation patterns. The high 
level of connectivity drives the network towards homogeneity. The 
agent accumulated locally will quickly be transported to other loca- 
tions, effectively reducing the local activation and preventing the 
nodes from activating. Only when the external feedback is strong 
enough to counteract this dissipation can the network support sus- 
tained activation (Fig. 4), leading to competition of the two processes. 
The network effectively supports only two modes of activity: random 
with formation of small domains (Fig. 4a locations ii and iii) and fully 
synchronous (Fig. 4a location i). 



We also investigated how the diffusion constant, D*, affects the 
degree of cooperation or competition in the system. To do this we 
subtracted the mean frequency, or mean phase coherence, for the 
active process case with no coupling (fi = 0.00) from that with strong 
coupling (P = 0.05), Fig. 5a, b. Comparing this baseline change for 
both metrics allows us to elucidate the extent to which the active 
process advances or delays changes in the frequency and MPC due to 
an increase in diffusion speed. We calculated the total area under the 
curves in Fig. 5a and b for all investigated values of y, Fig. 5c, d. 

The competition and cooperation patterns vary for both quantities 
and for both local and random network topologies. For local 
network connectivity, lower diffusion speeds lead to overall 
cooperation, whereas higher ones lead to competition. The general 
pattern of change in this case is similar for both frequency and MPC. 
The situation is different for networks with random connectivity. 
Here, the frequency shows competition for all of the speeds tested 
but its magnitude is non-monotonic with the diffusion speed 
(Fig. 5c). At the same time, for random topologies, the degree of 
cooperation and competition as reported by mean phase coherence 
changes largely monotonically as a function of the diffusion speed 
(Fig. 5d). 

The complex interaction of the two processes can be attributed to 
matched, for cooperation, or mismatched, for competition, patterns 
of agent spread, as supported by the two processes. The local spread 
of diffusion is defined by its speed, which determines the spatial 
range and the time scale within which the local agent level is affected. 
Network connectivity, on the other hand, determines the ratio of the 
local to global signal spread through the system. While the transmis- 
sion speed here is largely instantaneous, the effective speed of signal 
spread depends on actual activation levels of the nodes; that is, i.e. 
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Figure 3 | Formation of global activity patterns in networks with local connectivity, (a) Changes of the mean phase coherence as a function of active and 
passive transport amplitudes, (b) Mean phase coherence as function of active transport amplitude for y = 0.10. (c) Mean cross-correlation between pairs 
of nodes as a function of distance within the network, (d) Examples of network activation at a time t (the difference between the current time and the time 
of the last nodal activation), and histograms of pairwise cross-correlation values for different active transport amplitudes (as denoted in a). 



single activation of a node by the active process usually is not enough processes cooperatively may lead to activation of subsequent net- 
to trigger activation of a given site and is contingent on further agent work sites far from the original one. This however, due to the 
increase from other sources. Thus if the spread of passive and active intrinsic properties of diffusive signal spread, may only happen for 
processes match in their spatial scope and the timescale, then the two local network topologies. 
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Figure 4 | Formation of global activity patterns in networks with random connectivity. The individual panels correspond to those in Fig. 3. 
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Figure 5 | Cooperation and competition of active and passive processes as a function of changes in frequency and mean phase coherence (MPC) as 
functions of diffusion speed, D*. (a) Change in frequency between [5 = 0.05 and the baseline (/? = 0) as a function of passive transport strength, (b) 
Change in MPC between /? = 0.05 and the baseline (/? = 0) as a function of the passive transport strength. In panels a and b positive change denotes 
cooperation, while negative values denote competition, (c) Total increase and/or decrease in the frequency across all values of y due to the active process 
(calculated as a area under the curves in (a)), (d) Total increase and/or decrease in the MPC across all values of y due to the active process (calculated as a 
area under the curves in (b)). 



Competition arises when the activation of the nodes by one pro- 
cess is impeded by their earlier activation by the other process. This is 
predominantly (but not only) observed for random topologies. In 
this case the spatially random activation of nodes can be driven by 
two processes - noise and active spread from other sites through 
random connectivity. While both processes are independent of each 
other, activation of a given site (due to noise or active spread) and 
subsequent local activation of its neighbors due to passive agent 
spread through diffusion precludes instantaneous activation of that 
site due to the nodal refractory time since a node will not fire within a 
predefined time window after its activation. Thus the two processes 
compete for access to available nodes limiting the overall activation 
frequency. 

Discussion 

We have shown that network topology together with feedback from 
the external environment in which the network is embedded, may 
constitute a determining factor for the creation of sustained large- 
scale spatio-temporal patterns in a excitable system. We showed that 
local network topologies may cooperate with passive environmental 
transport to create globally propagating excitation through the sys- 
tem, whereas for random network topologies the two processes com- 
pete, increasing significantly the threshold in magnitude of the 
environmental transport needed for global patterns to occur. These 
results should be of importance in any system where the interaction 
of the network with its embedding environment can not be neglected. 
The work was supported by NSF CMMI grant no. 1029388. 

Methods 

Model. The spread of the agent is governed by two mechanisms: active spread 
through network connectivity, and passive environmental diffusion. The model 
describing the network nodes is very simple and is limited to charging and 
discharging elements that communicate through rapid equalization of the 
concentrations of the active agent. Such elements are relatively common in biology, 
and the model presented here roughly corresponds (to but is not limited) to astrocytic 
interactions. Passive diffusion is a natural choice for interactions through the 
environment as this is a common mode of environmental transport. Selecting 
diffusion as the interaction mechanism in the network requires the addition of a 



release like mechanism, which is accomplished by having an agent flush which is 
assumed to originate from an internal storage (Fig. 1). 

The equation describing the nodal dynamics is a modified integrate-and-fire 
model 14 : 



dXi 
~dt 



- oX f + I^eM + IflushM + (Xj - X, ) + y ]T Dy (r,t- Tj) ( 1 ) 



Xi is the amount of active agent within the z'-th node. The element time constant is set 
to z E = 3, while the nodal leak constant is set to a = 1. The I noise j is an instantaneous 
sub-threshold release defined as 



Incise At) =J2 I ^ H ( t ~ ~ H ( f ~ ~ AT )] • 



(2) 



Here, His the Heaviside function, the noise amplitude l n = 2.1; t nji denotes a specific 
noise instance; AT = 1 is the pulse duration. The noise arrives randomly at a given 
node with probability p„ = 0.01. Ijj us h,i describes the internal auto-release current 
when the level of the agent reaches the threshold, X { = 1 : 



Iflush,i(t) =Iat 



(3) 



where x { denotes the timing of the last threshold crossing in simulation steps. The size 
of the internal release is set to I amp = 1.5. The wave form is modeled as the difference 
between a slow, Q = 300, and a fast, C 2 = 30, exponential. Note that Ifl ush has a strong 
positive phase shortly after its activation which decays away over the duration of its 
refractory period. The refractory period, 1000 timesteps, prevents continuous 
activation of the system. 

The final term on the RHS describes the effect of the external, or environmental, 
excitation arriving at each network element, y is the coupling amplitude, and Dy(r, t — 
Tj) is the solution of the two dimensional diffusion equation at a distance r and time t 
in the space embedding the network. 



Dij{r,t 



N 



' 4nD*(t- 



-H)e-C(^) 



(4) 



For most of the simulations (except figure 5) N = 100 is the amplitude of neuro- 
transmitter release. The diffusion constant is D* = 0. 1, and the time decay constant of 
the agent in the external environment is £ = 0.01. For simplicity we assume that the 
element releases the agent into the physical space only when its value is supra- 
threshold. 

We use a Small World paradigm 15 to vary the network connectivity; initially the 
connections are set locally within given radius and then randomly rewired with 
probability p. The nodes are coupled through diffusive coupling with efficacy fi, with 
the connection from element; to element i denoted by Ay- taking the value 0 or 1. This 
paradigm allows for the range of network connectivity structures, from local to 
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random connections, to be controlled by the rewiring parameter, p. The strength of 
the active transport (f$ = 0.00-0.08) is limited by the amount of excitatory agent that 
it is capable of sending out at any moment, as well as how much a given element 
retains information of its past. The overall dynamics of the system strongly depends 
on the amount of the diffusive agent and its feedback interaction with the nodes. The 
parameters are chosen so that small network and environmental interactions result in 
a quiescent network while large environmental feedback results in hyperactivity of the 
nodes. The parameters near this transition were of interest (y = 0.08-0.14). 

Simulations. The model was simulated on a two dimensional network consisting of 
1600 oscillators on a 40 X 40 grid with periodic boundary conditions. The active 
connections on the network are initially assigned a radius R = 2. The network was 
initialized with no activity or charge present and allowed to evolve for 50,000 
timesteps using Euler's method. The initial 10,000 timesteps were removed from the 
analysis to allow the network to transition away from its (random) initial conditions. 
The figures displayed are analyzed over nine runs of the simulations for each set of 
parameters (ft, y,p). 

Analysis. The mean nodal frequency is calculated as the number of times any 
oscillator surpasses the charging threshold divided by the product of the number of 
oscillators and the timesteps in the simulation and then normalized. A frequency of 
one corresponds to a maximal raw firing frequency which is limited by the refractory 

time, j ' m a X ~ ^qqq • The units are not tied to actual timescales due to the generality of 

the model. 

The mean phase coherence 16 measures the degree of locking between nodal acti- 
vations, and is calculated from the relative nodal activation times. The pairwise mean 
phase coherence is defined by: 



1 N 



(5) 
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Where is the instantaneous phase of the j-th node relative to the z'-th node at fs 
activation at time, t^ k , 



(l>i(tj,k) =2n 



U,k,2 — tj,k 
U,k,2 — U,k,l 



(6) 



Here t^ k is the A;-th activation time of node j and t^i, U,k,2 are bracketing activation 
times of z'-th node such that t iykA — %k — U,k,i- The network wide mean phase 
coherence is the arithmetic average of the pairwise mean phase coherence. 
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